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ABSTRACT 

In the far future of an accelerating ACDM cosmology, the cosmic web of large- 
scale structure consists of a set of increasingly isolated halos in dynamical equilibrium. 
We examine the approach of coUisionless dark matter to hydrostatic equilibrium using 
a large N-body simulation evolved to scale factor a — 100, well beyond the vacuum- 
matter equality epoch, Ooq = 0.75, and 53h~^ Gyr into the future for a concordance 
model universe (fim = 0.3, f^A — 0.7). The radial phase-space structure of halos — 
characterized at a ^ acq by a pair of zero- velocity surfaces that bracket a dynamically 
active accretion region — simplifies at a ^ lOaoq when these surfaces merge to create 
a single zero- velocity surface, clearly defining the halo outer boundary, rhaio, and its 
enclosed mass, Mhaio- This boundary approaches a fixed physical size encompassing 
a mean interior density ~ 5 times the critical density, similar to the turnaround value 
in a classical Einstein-deSitter model. We relate Mhaio to other scales currently used 
to define halo mass (M200, -^vir, -^isofe) and find that M200 is approximately half of 
the total asymptotic cluster mass, while Migob follows the evolution of the inner zero 
velocity surface for a £ 2 but becomes much larger than the total bound mass for 
a ^ 3. The radial density profile of all bound halo material is well fit by a truncated 
Hernquist profile. An NFW profile provides a somewhat better fit interior to r2oo but 
is much too shallow in the range r2oo < r < rhaio- 
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1 INTRODUCTION 

In cosmologies in which density fluctuations are seeded by an 
early inflationary period (e.g., Kolb & Turner 1990), gravity 
acts during the era of cold dark matter domination to evolve 
a "cosmic web" of non-linear structure (e.g.. Bond, Kofman, 
& Pogosyian 1996) that can be approximately described as a 
set of roughly spherical halos, each characterized by a mass 
M. Within a given cosmology, the spatial density, clustering 
properties, and internal structure of such halos evolve with 
time in ways that are becoming increasingly well understood 
(e.g., Cooray & Sheth 2002, Kravtsov et al 2004). 

Oddly, the central element of this picture, the halo mass, 
remains poorly deflned. During the growth phase of the web, 
the mass distribution of a halo smoothly connects to the 
cosmological background of adjoining filaments, sheets, and 
voids. Any given halo is a mix of 'old' halo material that may 
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be near hydrostatic and virial equilibrium and 'new' mate- 
rial gained from recent accretion. Although a radial gradi- 
ent in the ratio of these two components is present, no clear 
edge that would allow a unique definition of mass separates 
them. In contrast, analytic collapse models based on spher- 
ical symmetry possess a well-defined outer shock or caustic 
surface (Bertschinger 1985; Filmore & Goldreich 1984) that 
emerges at a characteristic density. These models have moti- 
vated several mass definitions based on the enclosed density. 
White (2001) offers a recent discussion of these and other 
mass definitions applied to large-scale structure simulations. 

In contrast to the present epoch, the far future of a 
ACDM universe provides an opportunity to cleanly define 
halo mass. In the relatively near cosmological future, merger 
activity comes to an effective end and dark matter halos 
evolve toward a dynamically quiet state (Busha et al. 2003, 
hereafter paper I; Nagamine & Loeb 2003; see also the re- 
views of Adams & Laughlin 1997; Cirkovic 2003). In this 
paper, we show that the phase-space configuration of dark 
matter halos reaches a well-behaved asymptotic state char- 
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acterized by a single zero-velocity surface that uniquely de- 
fines the halo edge. Essentially all material internal to this 
surface is bound to and equilibrated within the halo, whereas 
material outside this surface is expanding away with the lo- 
cally perturbed Hubble flow. In §2, we describe the simu- 
lation used and the mass measures employed in our analy- 
sis. Results are given in §3, followed by brief summary and 
discussion section (§4). Our mass and radius measurements 
assume a Hubble constant Hq = 100ft km s~^ Mpc~^ with 
h = 0.7. 



2 SIMULATIONS AND MASS MEASURES 

We present results based on a simulation run with L- 
Gadget, a specialized version of the N-body code GAD- 
GET (Springel, Yoshida, & White 2001). L-Gadget simu- 
lates only coUisionless matter in a manner that optimizes 
cpu and memory resources. Our simulation models a patch 
of flat space in a ACDM universe with current matter den- 
sity ^Im = 0.3, vacuum density Qa = 0.7, and power spec- 
trum normalization erg — 0.9, values consistent with current 
observational data (Spergel et al. 2003). 

The dark matter in a periodic cube of side length 
200/i~^Mpc is modeled by 256^ particles of mass 3.97 x 
lO^°/i"^M0. The simulation was started at redshift z = 19 
(scale factor a — 0.05), evolved through the present epoch 
(a = 1), and continued forward to a = 100. A gravita- 
tional softening parameter fixed at 40/!.~^kpc (in physical 
units) was used throughout the computation. The simula- 
tion was run on 8 dual-cpu nodes of a Beowulf cluster at 
the University of Michigan. A total of 300 outputs equally 
spaced in log(a) were stored from the run. As the scale fac- 
tor increases, the universe evolves from being matter dom- 
inated to vacuum energy dominated. We define 

f-^cq SI'S the 

epoch at which the energy densities in the two components 
are equal. For the chosen model, this transition occurs at 
acq = (n™/nA)'/^ = 0.75. 

We identify dark matter halos using a standard per- 
colation (or friends-of-friends, EOF) algorithm with linking 
length 0.15 times the inter-particle spacing. The halo cen- 
ter is identified as the most bound particle of the result- 
ing group, and the halo velocity is defined as the center of 
mass velocity of the linked set of particles. At the end of the 
simulation this algorithm identified about 2900 halos with 
at least 500 particles. The largest halo at that time con- 
tains 83,600 particles, equivalent to 3.3 x W^^h~^MQ. In 
the analysis below, we measure ensemble properties using 
the 400 most massive halos. 

The FOE algorithm defines halo centers about which we 
measure several enclosed mass scales. The first mass scale 
is M200, the mass contained within a sphere of radius r2oo 
that encloses a mean density 200 times the critical density 
Pc(a) at the epoch of interest. In addition, a so-called virial 
mass Mvir (within r^ir) is defined to enclose a mean density 
Ac (a) times the critical density, where Ac (a) is an epoch- 
dependent threshold based on a simple, spherical collapse 
model (e.g.. Eke, Cole, & Frenk 1996). Ac reaches a con- 
stant value of 20.4 as flmia) — » 0. A third mass scale Mjsob 
is defined by a mean enclosed density of 180 times the back- 
ground mass density pm{a). 



3 RESULTS 

Previous studies (paper I; Nagamine & Loeb 2003) have 
shown that large-scale structure quickly approaches a sta- 
ble configuration in the future deSitter phase of a ACDM 
cosmology. Mergers and accretion slow dramatically after 
the scale factor exceeds a = 2 — 3, signifying the end of 
the dynamically active stage of the cosmic web. The change 
in dynamical activity, from active to absent, is apparent in 
the reduced phase-space density f{r,Vr) of halos as shown 
in Figure Q This figure plots the proper radial velocity 
Vr against distance from the halo center r for an ensem- 
ble average of the 400 most massive halos using r2oo and 
V200 = a/ GM200 /r200 as scale variables. 

Figure ^ shows the conditional probability p{v,. \ r) = 
f{r,Vr)/p{r) for these 400 halos at a = 1. The solid line 
gives the mean radial velocity, while the grey scale regions 
delimit velocities containing 40%, 60%, 80%, 95%, and 99% 
of the material at each radius. This ensemble average phase 
space density can be divided into three principal regions: i) 
an inner hydrostatic core ((wr) = 0) that is relatively well 
relaxed; ii) an intermediate accretion envelope {{vr) < 0) 
containing two opposing streams of material, one on its first 
inward journey toward the halo center and the other passing 
outward after pericentric passage; and iii) an outfiow region 
{(vr) > 0) dominated by the locally perturbed Hubble fiow. 
Nearly all of the material in regions i) and ii) is gravitation- 
ally bound to the halo (a modest fraction of material within 
r2oo can be scattered out of a halo, see Fig. 3 of paper I) 
while essentially all of region iii) is unbound. The character- 
istic scale that separates the outfiow and accretion regions is 
often called the turnaround radius rta (Gunn & Gott 1972), 
a term motivated by spherical models of expanding mass 
shells. We call the radius that separates the two interior 
zones the hydrostatic radius rhs and measure its value using 
a threshold condition for the binned, mean radial velocity. 
Starting from the interior, we identify rhs as the minimum 
radius at which |(fr)/'y20o| > 0.1. The turnaround radius is 
identified in the same manner, identifying the maximum ra- 
dius at which |(i'r)/v20o| < 0.1. At a = 1, the values of these 
characteristic radii are rhs = 0.70r2oo and rta = 3.3r2oo and 
are marked by the vertical dashed lines in Figure^. 

The averaging process used to create the density in Fig- 
ure^ blends the effects of substructures within individual 
halos, especially in the hydrostatic region. At a given radius, 
this region has a nearly Gaussian distribution of radial ve- 
locities with zero mean, signatures of hydrostatic and virial 
equilibrium. Due to the presence of pairs of massive clus- 
ters, the outfiow region for the ensemble profile has a rather 
broad dispersion. The fact that the distribution of radial ve- 
locities for the outflow region is more sharply peaked than in 
the virialized region indicates that the surrounding regions 
generally contain halos that are less massive than those of 
our high-mass selected sample. 

The situation in the far future is markedly different from 
the present. Eigure^D shows the ensemble phase-space den- 
sity at a = 100. The most striking changes are the disap- 
pearance of the accretion region and the dramatic cooling of 
the outflow region. Here, rhs and rta have merged to form a 
single zero- velocity surface at rhaio ~ 4.5r2oo, represented by 
the dashed vertical line. In addition, relatively few nearby 
halos are present to disrupt the outflow stream. Neighboring 
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Figure 1. The distribution of dark matter radial velocities as 
a function of distance from the halo center. The top two panels 
show the conditional phase-space density p{vr\r) as a function of 
radius for the ensemble of 400 largest halos at the present epoch 
(a) and for the future when a = 100 (b). The solid line shows 
the mean velocity as a function of radius; the grey scale indicates 
the regions enclosing 40, 60, 80, 95, and 99 percent of the particle 
population as specified by p{vr\r); the vertical lines represent the 
zero-velocity surfaces. Panel (c) shows the mean radial velocity 
for an ensemble of halos at epochs a = 0.34, 0.59, 1.0, 1.8, and 
4.8, with the bold line representing the function at a = 100. 

halos that existed at a = 1 have either merged or been pulled 
away in the deSitter expansion. Although not shown here, 
the phase-space density of the single most massive halo at 
a = 100 is nearly identical to the ensemble averaged profile. 
Because mergers cease long before a = 100, the substruc- 
ture present within any individual halo has several dynam- 
ical times to relax via tidal stripping, phase mixing, and 
dynamical friction (see paper I for a comparison with the 
phase-space distributions for individual halos). 

The elimination of the accretion region at a = 100 is 
accompanied by an expansion of the hydrostatic region. Fig- 
ure^ shows the evolution of the average i'r/v200- The bold 
curve shows the final profile at a = 100; the thinner curves 
depict the epochs a = 0.34, 0.59, 1.0, 1.8, and 4.8. The ex- 



tent of the accretion region decreases with time as rhs/r2oo 
grows and rta/r'2oo slightly shrinks. Note also that there is 
an expulsion epoch at a = 4.8, shortly after mergers have 
ended. At this time, Vr/v2oo is somewhat larger than its 
asymptotic value in the range r/r2oo = 3 — 8 because the 
unbound particles from the final mergers are being ejected. 

The time evolution of this transition is presented in Fig- 
ure|21 As before, we show ensemble behavior of the 400 most 
massive halos selected at a = 100, with halos at earlier epochs 
restricted to the most massive progenitors of this final pop- 
ulation (where most massive progenitor is defined to be the 
halo with the largest Mfof that contributes at least 20% 
of its particles to its descendant). For each halo and epoch, 
we identify the physical values of the characteristic radii, 
(rhs, rta, 7'2oo, fvir, Tisob) aloug with the respective enclosed 
masses. To determine rhs and rta for each halo, we first time 
smooth the individual profile by co-adding the halo config- 
uration over seven consecutive outputs. We then measure 
rta using a linear extrapolation of the outflow over a factor 
of six in Vr, starting at the radius r = lOMpc and working 
inward. Applying a threshold (as above) instead of extrap- 
olating produces somewhat larger values for rta, but causes 
only a small change in the values for Mta. We choose the ex- 
trapolation method because it provides less noisy estimates 
during the early phase of active halo growth. The hydro- 
static region, r^s is measured using a threshold technique on 
the mean radial velocity measured in radial bins, identical to 
the method used for Figure^ We use logarithmically spaced 
bins, 30 per decade, an d identify rhs as the radius at which 

\Vr\/vta > 0.1, Vta = y/ GMta/rta. Typically, Vta ~ 0.7U200- 

Figure 121 shows the time evolution of the mean physical 
radii and enclosed masses for the different measures. During 
the early Einstein-deSitter phase, a < acq, the sizes defined 
by mean interior densities r2oo, rvir and rigob are similar, 
whereas the hydrostatic boundary rhs hes somewhat interior 
to r2oo and rvir. The turnaround radius rta hes well beyond 
these scales while the growth of linear perturbations remains 
robust. 

The interior scales diverge at epochs a > a^q when 
vacuum effects dominate. As linear growth stagnates, large- 
scale mass density fluctuations 5p/p become frozen into the 
background. Continuity of the density field implies that the 
mean size rigot. (defined relative to the mean matter density) 
becomes constant in the comoving frame, so the physical 
size grows exponentially in time. In contrast, the ensemble 
mean r2oo rapidly approaches its asymptotic value, follow- 
ing r2oo(a) = r2oo,oo [1 — exp(— (a/aoq)^'^)] at late times, 
with r2oo,oo = 1.2 Mpc for this high-mass sample. This fit is 
shown for a > Oeq by the dot-long-dashed line in Figure 

The value of rvir relaxes less quickly, due to the decline 
in the variable threshold Ac (a) as the scale factor increases. 
The enclosed masses M200, Mvir and Migob grow monotoni- 
cally until reaching 99% of their asymptotic limits at a = 2.9, 
8.6 and 13.6, respectively, corresponding to ages of 32, 52 
and 60 Gyr. The innermost mass scale A'hoo experiences a 
slight decline at late epochs, oc (0.00045 ± 0.00015)ln (a), but 
it is not clear whether this drift is physical or numerical. 

The hydrostatic boundary r^s tracks rigot. until a ~ 2, 
then slows its growth and relaxes toward an asymptotic 
value. The turnaround radius also grows rapidly until a ~ 
1.5, then declines until a ~ 5, and increases slowly there- 
after. The decline in size and enclosed mass from the peak 
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Figure 2. Comparison of mass and radial scales for the progen- 
itors of the 400 most massive halos identified at the end of the 
computation. The four panels show the scale factor dependence 
of (a) mean masses; (b) mean physical sizes; (c) enclosed den- 
sities relative to the critical value. The line styles indicate the 
different mass measures: Mta (bold); Afj^g (dashed); M200 (solid); 
Mvir (dotted); and Migot (dot-dashed). The line styles for the ra- 
dial scales are analogous. In panel (b), the dot-long-dashed curves 
show the asymptotic form of r2oo (see text) for a > acq = 0.75. 



until a ^ 5 arises from a change from infall to outflow in 
the accretion regime lying between rhs and rta, as indicated 
by the outflow enhancement at this epoch in Figure At 
o~ 1, accretion just outside the hydrostatic boundary drops 
dramatically. After a crossing time, the modest mass frac- 
tion of this accreted material that is scattered to unbound 
orbits emerges at radii r > rhs with positive radial velocity 
at a ^ 2 — 5. This material escapes from the halos, but we 
have not investigated its ultimate fate. Presumably, some 
fraction emerges as the cores of stripped sub-halos that re- 
main bound, whereas the remainder may emerge sufficiently 
tenuous and hot that it never collapses into a halo of cosmo- 
logically interesting mass. Future studies using higher reso- 
lution experiments are needed to address this question. 



For test particles in Hubble flow around a halo of mass 
M, the outer zero- velocity surface location can be estimated 
using a Newtonian binding energy, with the result (eq. [11] 
of paper 1) r ~ (A//lO^^M0)^''^Mpc. Applying this estimate 
using the mean asymptotic halo mass in Figure yields a 
value of 5Mpc, close to the asymptotic mean size shown in 
Figure|5j3. We caution that the values of rhs and rta are sen- 
sitive at the ~ 10% level to the choice of threshold and/or 
the interpolation scheme used to locate the zero-crossing. 
However, the steep behavior of the density profile near the 
boundary (pocr~^, 7^4) leads to a much smaller uncer- 
tainty ( 2%) in enclosed mass. 

The merging of the turnaround and hydrostatic scales 
signals the end of high mass structure formation, and there- 
fore the end of the growth phase of the cosmic web. There- 
after, halos maintain a fixed physical size and the morphol- 
ogy of the cosmic web, when viewed in the comoving frame 
at a density threshold near critical, becomes less of a web 
and more an increasingly fine spray of droplets (see paper 
I and Nagamine & Loeb 2003). From Figure the ratio 
Mhs/A/ta reaches 99% of its asymptotic value at a = 7.4 
(a/aoq = 10). At epochs beyond lOttcq, we can clearly define 
the (ultimate) halo mass as Mhaio = M^s ~ Mta, i.e., the 
mass enclosed within the single zero-velocity surface is the 
ultimate halo mass. 

Figure shows that the turnaround/ultimate halo 
mass is defined by a density threshold phaio relative to crit- 
ical that lies at all times within a factor of two range span- 
ning ~(5 — lO)pc(a). The threshold at early times lies close 
to the canonical 97r^/16 = 5.5 value expected from spheri- 
cal collapse in an Einstein-deSitter universe (Peebles 1980). 
After Ocq, the ratio Phaio/pc(a) first drops while the mean 
halo mass peaks, then climbs to 10 during the period 
a /acq = 2 — 10 when the mean halo mass drops. At late 
times, phaio is declining weakly, and its asymptotic limit, al- 
though not well determined by this simulation, is only about 
12% shy of the canonical value. 

While the gap between rigob and rhaio expands exponen- 
tially in the far future, the enclosed mass ratio converges to 
Afisob/A'^haio ~ 1.35. While thresholding with respect to the 
background mass density extends beyond the halo edge at 
late times, the scales defined relative to the critical density 
pick out mass shells interior to the halo: Mvir = 0.89Mhaio 
and Af2oo = 0.52Afhaio. These particular values are sensitive 
to the asymptotic form of the radial density distribution. 

Figure shows the mean radial profile obtained at a = 
100 from the stacked ensemble of the 400 most massive halos. 
We use r2oo as the scale radius, but results are similar when 
other characteristic scales are used. The light solid line shows 
the mean profile while the heavy solid line shows the profile 
for material bound to each halo, using specific energy E/m — 
w^/2 - GM{< r) /r + (47r/3)GpAr^ (with v and r the proper 
velocity and radius, respectively). The bound material is 
well fit over the entire halo volume by a truncated Hernquist 
(1990) profile 
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with characteristic radius rc = 0.62r2oo and truncation scale 
rhaio = 4.6r2oo. The latter measure of halo size agrees ex- 
tremely well with the value of 4.5r2oo obtained from the 
ensemble average zero- velocity surface. This fit has a least- 
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Figure 3. The asymptotic form of the density distribution for 
dark matter halos. The solid curves show the density profile for an 
ensemble average of the 400 most massive halos in the simulation 
at scale factor a = 100. The upper solid curve shows the total 
density as a function of radius, whereas the lower solid curve 
includes only the bound particles. The dotted curve shows the 
best fit NFW profile. The dot-dashed curve shows the best fit 
Hernquist profile, and the dashed curve represents a truncated 
Hernquist profile (see eq. [Til. 

square error of about 12%. For comparison, an NFW profile 
p(r) = 4ps(r/rs)"^[l + (r/r^)]"^ (Navarro, Frank, & White 
1997), obtained by fitting within r2oo is shown in Figurej^by 
the dotted line. The profile with rs =0.25r2oo fits the inner 
regions of the halo well, but it is much too shallow beyond 
r2oo, overestimating the halo density at rhaio by an order of 
magnitude. Note that simple scaling for the evolution of the 
concentration proposed by Bullock et al (2001), Cvir{a,) ~ a, 
breaks down as the growth rate becomes constant and the 
profile out to rvir becomes fixed after a few Hubble times. 



4 SUMMARY 

From a large N-body simulation that follows the long term 
evolution of coUisionless structure in a ACDM cosmology, 
we have examined the ultimate approach to equilibrium 
for a sample of 400 massive halos with final mean mass 
3.5 X 1O^*M0. During the matter-dominated era (a aeq = 
0.75), the radial phase-space structure of the halos is com- 
plex, consisting of an inner hydrostatic region, an inter- 
mediate accretion zone, and an exterior region that ex- 
pands with the perturbed Hubble flow. During the interval 
Oeq < a < lOaeq, accretiou shuts down and the intermedi- 
ate region is briefly dominated by outflow rather than infall. 
For a ^ iOacq, the intermediate region disappears, and the 
hydrostatic and turnaround scales merge to form a single 
zero-velocity surface that provides an unambiguous defini- 
tion of halo mass. 

The existence of multiple mass scales commonly used 
in the literature to describe clusters is a direct reflection of 
the complexity of the accretion region during the growth 



phase of structure. This work illuminates the evolving rela- 
tionships between the hydrostatic and turnaround scales of 
halos and scales defined by mean interior density thresholds. 
Although thresholding with respect to the background den- 
sity (Misob) has advantages for calculating the halo space 
density (Jenkins et al. 2001), its use to describe future struc- 
ture is compromised by our finding that Misob exceeds the 
asymptotic halo mass for a ^ 2. Masses defined by thresh- 
olds relative to the critical density converge to well-defined 
values. For a ^ lOttcq, we find A/haio = 1.1 Mvir = 1.9 A/200. 

Using only bound material, the ensemble mean density 
profile of the 400 most massive halos is well fit by a modified 
Hernquist model, with scale radius 0.62r2oo and truncation 
radius 4.6r2oo (eq. 0). The origin of this particular form, 
as well as other issues such as its dependence on halo mass 
and its extension to an ellipsoidal description, remains to be 
investigated. 
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